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ABSTRACT 


Observing gravitationally lensed objects in the time domain is difficult, and well-observed time- 
varying sources are rare. Lensed gamma-ray bursts (GRBs) offer improved timing precision to this 
class of objects complementing observations of quasars and supernovae. The rate of lensed GRBs is 
highly uncertain, approximately 1 in 1000. The Gamma-ray Burst Monitor (GBM) onboard the Fermi 
Gamma-ray Space Telescope has observed more than 3000 GRBs making it an ideal instrument to 
uncover lensed bursts. Here we present observations of GRB 2108124 showing two emission episodes, 
separated by 33.3 s, and with flux ratio of about 4.5. An exhaustive temporal and spectral analysis 
shows that the two emission episodes have the same pulse and spectral shape, which poses challenges 
to GRB models. We report multiple lines of evidence for a gravitational lens origin. In particular, 
modeling the lightcurve using nested sampling we uncover strong evidence in favor of the lensing 
scenario. Assuming a point mass lens, the mass of the lensing object is about 1 million solar masses. 
High-resolution radio imaging is needed for future lens candidates to derive tighter constraints. 


Keywords: gamma-ray burst — gamma-ray transient source 


1. INTRODUCTION 


Strong gravitational lensing is a tool that serendip- 
itously enhances our observing capabilities and offers 
new opportunities to study the Universe (see e.g. Oguri 
2019). Gamma-ray bursts (GRBs) are energetic tran- 
sient sources at cosmological distances, involving rela- 
tivistic jets from stellar-mass black hole (BH) central 
engines. GRBs last from a fraction of a second to about 
1000 s and typically show non-thermal spectra (see e.g. 
Kumar & Zhang 2015; Beloborodov & Mészáros 2017, 
for reviews). Given that the distance scale of GRBs 
spans a wide range (up to redshift z < 9, Cucchiara 
et al. (2011)), a fraction of GRBs will show the imprints 
of strong gravitational lensing. 

Strong gravitational lensing produces multiple images 
of the same source. The images differ in their inten- 
sity, but importantly their spectral shapes will be the 
same. Similarly, in time-varying sources, the temporal 
profile will be the same but shifted in time for differ- 
ent images (Schneider et al. 1992). The temporal and 
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spectral invariance is the main defining feature of grav- 
itational lensing. Lensed images are separated by up 
to arcseconds which are clearly below the resolution of 
current gamma-ray detectors. Gamma-ray instruments, 
however, have excellent time resolution, and temporal 
structures can be recorded with unparalleled accuracy 
(Meegan et al. 2009). 

In one incarnation of the lensing scenario (Refsdal 
1964; Rodney et al. 2021) applied to GRBs, also called 
macrolensing, a single event triggers the same instru- 
ment twice, and can be separated anywhere from a few 
hours to decades. The two triggers will have similar 
lightcurve shapes and spectra. The delay between the 
events scales linearly with the lens mass. Lens candi- 
dates for this scenario have masses in the 105 — 101? Mo 
range. Objects in this mass range include supermassive 
black holes (up to few x10!? M), galaxies and clusters 
of galaxies. In practice, however, all traditional strong 
lensing time delay measurements are from galaxies or 
clusters of galaxies. GRB lensing rates are highly un- 
certain but somewhere on the order of 1 in 1000 GRBs 
should be affected (Mao 1992). After roughly 10,000 
observed GRBs during three decades, no convincing 
macrolensing candidate GRB pair has been found (Ne- 
miroff et al. 1994; Davidson et al. 2011; Veres et al. 2009; 
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Hurley et al. 2019; Ahlgren & Larsson 2020). The neg- 
ative result is likely a combination of two effects: First, 
GRB detectors typically in low Earth orbit will miss a 
sizeable fraction of GRB lens echoes (see, however Hur- 
ley et al. 2019; Hui & MoonBEAM Team 2021, for exist- 
ing and future all-sky instruments). Second, for weaker 
GRBs, with pulses close to the noise level, it is difficult 
to distinguish between the lensing scenario and just two 
unrelated but similar looking GRBs (Ahlgren & Larsson 
2020). 

In a different scenario, also called millilensing (Ne- 
miroff et al. 2001, because the expected separation be- 
tween the images is on the order of milli-arcseconds), 
the gravitational lens signature is imprinted upon the 
lightcurve of a single trigger. In this case, we have, e.g., 
two emission episodes with similar lightcurve patterns 
which can be separated by timescales spanning from a 
fraction of a second to a few minutes. 

Recently, there has been an increase in claims of 
millilensing events. Paynter et al. (2021) presented 
convincing evidence for lensing in the short duration 
BATSE GRB 950830. Mukherjee & Nemiroff (2021a) 
raised some concerns based on the inconsistent flux ra- 
tio between the two pulses. Yang et al. (2021) and 
Wang et al. (2021) independently argued that the likely 
short GBM GRB 200716C shows milli-lensing signa- 
tures. Kalantari et al. (2021) reported on a different 
GBM lensing candidate, the long duration GRB 090717, 
selected based on the analysis of the auto-correlation 
function. The claim of Kalantari et al. (2021) was chal- 
lenged by Mukherjee & Nemiroff (2021b) arguing that 
the two pulse shapes differ significantly. For the above- 
reported lensing candidates the flux of the first and the 
second emission episode is either at the same level or 
their ratio is S; 1.5. 

We present observations of the long duration 
GRB 2108124 and show that it is consistent with a grav- 
itational lensing scenario. It is the first lensing claim 
with a flux ratio = 3. We list multiple lines of evidence 
to support the lensing interpretation. We perform spec- 
tral and temporal analysis using Fermi-GBM data and 
complement it with additional data from INTEGRAL- 
SPI/ACS and Swift/BAT. 

In Section 2, we present the observations, followed by 
tests for gravitational lensing in Section 3. We discuss 
the power of the tests in Section 5 and conclude in Sec- 
tion 6. 


2. OBSERVATIONS 


Fermi-GBM triggered on GRB210812A (Figure 1, 
trigger number 650479626 / 210812699) on August 12°”, 
2021 at 16:47:01.014 UT (T0) (Fermi GBM Team 2021; 
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Figure 1. Background subtracted  lightcurve of 


GRB210812A. The data is summed over all GBM de- 
tectors with good coverage. The red lines show a subset of 
MCMC fits using the N1 pulse (see Equation 3 and Section 
4.1.1). 


Veres & Fermi-GBM Team 2021). Fermi-GBM consists 
of 12 Nal (referred to as n0, n1,..., n9, na and nb) and 2 
BGO (referred to as b0 and b1) detectors covering the 
entire unocculted sky in the 8-1000 keV and ~0.1-40 
MeV energy range, respectively. 

As reported by the automatic pipeline, the GRB lo- 
cation is RA= 39.7, dec= 69.7 degrees, with an error 
radius of 1.1 degrees (statistical only). At the trigger 
time, this position corresponds to a LAT boresight an- 
gle of 149 degrees. The GRB location was behind the 
spacecraft, meaning most GBM detector normals have 
a large angle to the source. Based on previous expe- 
rience (e.g. Connaughton et al. 2016), this type of ge- 
ometry results in a large number of detectors showing 
approximately equal count rates compared to ~3 detec- 
tors with dominating signal for typical GRBs with small 
LAT boresight angle. Upon visual inspection of the Nal 
detectors’ lightcurves in the 50-300 keV range, where 
GRBs are brightest, we find that all 12 Nal detectors 
detected the brighter first peak. Moreover, detectors 
nl and n6 through nb also detected the fainter, second 
pulse. Both of the pulses were detected by the 2 BGO 
detectors. 

Detectors n8, na, and nb showed the strongest signals. 
We used data from these detectors, along with bl, for 
spectral analysis. For the temporal analysis, we used de- 
tectors n1, n6 through nb, b0, and b1. For both spectral 
and temporal analysis, we used the 128 energy channel 
time-tagged event (TTE) data. We chose the pre-binned, 
8 energy channel (ctime) data to carry out the flux-ratio 
test. 
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The targeted search (Blackburn et al. 2015; Gold- 
stein et al. 2019) was designed to search for coherent 
sub-threshold signals (weak signals that did not trig- 
ger the instrument). As expected, we recovered both 
pulses of GRB210812A with high significance. The 
search also provides a location for the burst (RA= 40.5°, 
dec= 69.4°), consistent with the location of the auto- 
matic pipeline used for standard GRB analysis. We also 
find that the locations of the two pulses are consistent, 
meaning they do indeed belong to the same source. 


2.1. Other observations 


The location of GRB 210812A was outside the coded 
mask of Swift-BAT (Gehrels et al. 2004; Barthelmy et al. 
2005) at the time of the trigger. The first peak of 
GRB 210812A is clearly present in the continuous 4- 
channel data!,?, but the second peak is only discernible 
in the summed lightcurve. 

INTEGRAL/SPI-ACS (Winkler et al. 2003; von Kien- 
lin et al. 2003) detected GRB 2108124 and clearly shows 
the two-peak structure.? We calculate the time delay be- 
tween Fermi-GBM and ACS is due to the higher altitude 
of the INTEGRAL spacecraft to be dtacg = —0.396 ms. 
We applied this correction to the ACS lightcurve in Fig- 
ure 2. 


2.2. Temporal properties 


GRB 210812A consists of two pulses, separated by a 
quiescent period of about 30 s (Figure 1). Using the 
GBM targeted search, we found no emission between 
the two pulses. For background estimates, we fit a third 
degree polynomial based on quiescent segments of the 
lightcurve before, after and in-between the two pulses. 
We also report no discernible emission around 33 s after 
the second pulse, indicating that this is not a periodic 
source. We further note that low level emission dis- 
cernible in the summed lightcurve just before the start 
of the second pulse (T0+26.5 s to T0+29.6 s) is incon- 
sistent with coming from the location of GRB 210812A, 
and thus it is unrelated. 

GRB 210812A has a duration of Too = 39.9 + 3.6 s 
(10-1000 keV; Too marks the time interval between 5% 
and 9596 of the cumulative flux)*. Analyzed separately, 
the two pulses have consistent duration: the duration of 
the first pulse is 799, = 5.31 + 0.68 s while the second 


1 https:/ /www.swift.ac.uk/archive/reproc/00096363010/bat / 
rate/sw00096363010brtms.lc.gz 


? https:/ /www.swift.ac.uk/archive/reproc/00036688035/bat / 
rate/sw00036688035brtms.lc.gz 


3 http://isdc.unige.ch/-savchenk/spiacs-online/spiacs.pl 
4 https:/ /gcn.gsfc.nasa.gov/gcn3/30633.gcn3 
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Figure 2. GBM lightcurves in different energy ranges, 


summed over all detectors with good coverage. ACS and 
Swift lightcurves are also shown in the bottom two panels 
with the native resolution (black, 50 ms for ACS and 64 ms 
for BAT) and binned (red, 0.4 s for ACS and 0.512 s for 
BAT). 


pulse is 799,5. = 3.84 + 1.64 s long. Taking the first 
pulse as a separate GRB, we classify it based on the 
GBM Too distribution (Bhat et al. 2016; von Kienlin 
et al. 2020), as a likely long GRB originating from the 
collapse of a massive star, with a probability of 8796. 
Conversely, the likelihood of GRB 2108124 being a short 
GRB from a compact binary merger, based on the T99,1 
information, the observed distribution of Fermi-GBM 
GRBs and calculated consistently with Goldstein et al. 
(2017); Rouco Escorial et al. (2021) is 13 96. 
Autocorrelation function - We measure the time delay 
between the two pulses using the autocorrelation func- 
tion of the lightcurve between 10 and 1000 keV, with 64 
ms resolution. To accurately determine the autocorre- 
lation peak, we fit a 9th degree polynomial to the peak 
region. We estimate the uncertainty of the delay by 
adding Poisson noise to the original lightcurve (see e.g. 
Ukwatta et al. 2010; Hakkila et al. 2018, for a similar 
approach). After repeating the peak finding procedure, 
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Figure 3. Comparison of the lightcurves of the two pulses. 
The stronger first pulse is scaled down to match the second 
pulse. The fainter second pulse is shifted in time to match 
the first pulse. Dark and light regions are 1 and 2 sigma 
regions respectively. The figures represent: top-left: Nal 
(45-300 keV) + BGO (120-400 keV), top-right: NaI (22-800 
keV), bottom left: ACS (> 70 keV), bottom-right: BGO 
(120-400 keV). 


we take the uncertainty as the lo confidence region of 
the delays of the modified lightcurves and get: 


Ataor = 33.301017 s. (1) 


This is the first among the many time delay measure- 
ments between the pulses. We note here that this ac- 
curacy is typical of what one would expect for a lensed 
GRB separated by a longer timescale (macrolensing). 

The peak of the auto-correlation curve shows a 3.150 
excess over the smoothed curve calculated using the 
Savitzky-Golay filter. This excess satisfies the crite- 
ria of Paynter et al. (2021) for lensing, which classi- 
fies GRB210812A as a lensing candidate for further 
scrutiny. 

Spectral lag - The lag is a measure of the delay between 
high and low-energy photons (e.g. Norris et al. 2000). 
Typically, a GRB is harder initially, and the higher en- 
ergy photons (100-300 keV) arrive earlier than the low 
energy (25-50 keV) photons. This relation is also used to 
estimate the redshift based on the empirical correlation 
between lag and luminosity (Norris et al. 2000). For the 
first pulse we find that the lag is Tiag,1 = 221.6* 1224 ms, 
while the second pulse has: Tiag2 = —89.2733¢2 ms. 
The error on the second pulse is much larger, because of 
its weakness, especially in the 25-50 keV range. 


2.3. Spectral Analysis 


The spectrum of the first pulse is best fit by a power 
law with exponential cutoff, also known as the Comp- 
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Figure 4. Spectra of the two pulses. The second pulse is 
shifted to match the first. Shaded regions mark 1 and 20 
confidence regions. 


tonized model (see Table 1 for the parameters and Fig- 
ure 4 for the spectrum). The fluence in the 10-1000 keV 
range is F, = (80.1 + 0.3) x 1077 erg cm ?. 

'The second pulse is weaker, and it is fit equally well 
by a simple power-law and the Comptonized model. We 
chose the Comptonized model, because it is more phys- 
ical (the power law with photon index > —2 integrates 
to infinite energy). The difference in the goodness of fit 
measure (AC-stat) when going from the simpler power 
law (2 parameters) to the Comptonized model (3 param- 
eters) is ~6. This is just below the AC-stat ~ 8 used 
in e.g. Poolakkil et al. (2021) to select the more com- 
plex model. However, all the parameters of the Comp- 
tonized model are well constrained (Table 1), which jus- 
tifies the use of this model. The fluence of pulse 2 is 
Fy = (21.6 + 2.6) x 1077 erg cm-?. 

We can also compare the time-resolved spectrum of 
the two pulses. Because of the relative weakness of the 
second pulse, we chose to fit the simple power-law model 
in bins of 0.256 s. We show the temporal evolution of 
the power-law indices for the two pulses and a linear fit 
to both as a function of time (Figure 5). We shifted the 
second pulse by 33.3 s to highlight their similar behav- 
10r. 


3. INDICATORS OF LENSING ORIGIN 


Proving the gravitational lensing origin involves show- 
ing that the two pulses have identical pulse shapes, 
and their spectra are similar as well. It is conceivable 
that some physical mechanism produces two pulses with 
identical pulse shapes and spectra with no lensing in- 
volved. Without imaging the sources based solely on 
the gamma-ray observations, a lensing scenario is diffi- 
cult to prove beyond the shadow of a doubt. 
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Interval Epeak photon index C-stat/dof photon flux Energy flux 
(s) (keV) (ph cm~? s~') (1077 erg cm ? s^) 
-0.256 to 4.096 | 324+28 —0.88 +0.08 373.3/310 10.35 + 0.46 18.4+0.7 
33.024 to 37.376 | 283 +90 —1.10 +0.24 295.3/304 3.48 + 0.49 4.96 + 0.59 


Table 1. Spectrum of the two pulses fit by a power law with exponential cutoff. 
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Figure 5. Temporal evolution of the spectral power law 
index compared for the two pulses. The slopes of the fit- 
ted linear functions are indicated in the legend and shaded 
regions mark 1c confidence intervals. 


In this section, we will show, however, that the 
gamma-ray properties of GRB 2108124 pass all the tests 
in the literature for a lensing origin. Furthermore, we 
apply the Bayesian evidence criterion that can select 
among models, and this method indicates strong evi- 
dence in favor of the lensing scenario. 

The basic idea for establishing the statistical likeli- 
hood of lensing in the temporal domain is comparing 
two scenarios. First, we assume no lensing. We fit the 
pulse model and allow every parameter to vary freely. In 
the alternative scenario, where lensing is assumed, the 
second pulse is forced to have the same shape as the first 
one and differ only by a normalization factor and a time 
delay. 

We use two types of pulses (Norris et al. 1996, 2005). 
The first pulse with 5 parameters, referred to as “N1”, 
is defined as: 


— (fmax =t)” " 
e or itt 

Iui(t|A, tax, Or; Ca; V) =A - (Spm )" add 
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Here A is the amplitude at peak time, tmax. c, and 
cq mark the rise and decay timescales of the pulse and 
v is a shape parameter. 


The second pulse with 4 parameters, ^N2", has the 
following temporal dependence: 


t—A 


Ina(tlA, £, A, 7) =A ot = 7 (2) 


where A is the amplitude, € is the asymmetry parameter, 
A is the start time of the pulse, and 7 is a duration 
parameter. 


3.1. Indirect evidence 


We explore a few properties of GRB 210812A that are 
not direct proofs of lensing, but they are necessary to 
any such claim. 

Hard-to-soft evolution: Gamma-ray bursts typically 
become softer with time (e.g. Kaneko et al. 2006). 
This trend can be observed in the individual pulses of 
GRB 210812A as well (see Figure 5). For GRBs in gen- 
eral, the hard-to-soft evolution can be observed even 
across pulses. Specifically, for GRBs that show two 
pulses separated by a quiescent period, the second pulse 
is, in general, softer (e.g. Lan et al. 2018; Zhang et al. 
2012). We find that the second pulse of GRB 210812A 
has the same spectral shape within errors or, subse- 
quently, the same hardness, which is uncommon for typ- 
ical GRBs. 

Leading pulse is brighter: In case of simple lens mod- 
els, the light ray traveling closer to the lens arrives later. 
It has a lower magnification than the first arriving light 
ray with the larger impact parameter (Krauss & Small 
1991). We find that the first pulse is indeed visibly 
brighter and thus consistent with the expectation from a 
lensed source. We note that this criterion may not hold 
for complex lens models (Keeton & Moustakas 2009). 


3.2. Spectrum of the pulses 


For the spectral analysis we first selected the interval 
containing the second pulse by visually identifying con- 
tiguous temporal bins with significant signal. For the 
first pulse, we selected a source interval which is the 
same length as the second pulse (see Table 1 and Figure 
2). The spectral fits of the two pulses yield consistent 
spectral shapes within errors: the peak of the energy- 
per-decade or vF, spectrum is Epeak,ı = 324 + 28 keV 
and Epeak,1 = 283 + 90 keV. To compare the two spec- 
tra, we plot the 1 and 2 o confidence regions of the 
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spectral shapes (Figure 4), accounting for the correla- 
tions between the parameters. We multiply the second 
pulse by the fiducial 3.5 number to show that the two 
spectral shapes are consistent. 


3.3. Count ratio test 


If the pulses are gravitationally lensed, the ratio be- 
tween pulses should not depend on energy. Mukherjee 
& Nemiroff (2021a) investigated the count ratio of the 
two pulses as a function of energy for GRB 950830 and 
found a = 2e inconsistency. This test has the advantage 
that it is independent of the particular spectral shape. 

For GBM we used the 8-channel ctime data to carry 
out this test on GRB 210812A. We considered all the de- 
tectors where the second pulse was visible in any channel 
(nl, n6 through nb, energy channels 1 through 6, and 
b0 and bl, energy channels 0 and 1), and data from 
ACS and BAT (25-350 keV). We find that the ratio 
of the pulses in all the channels and all three instru- 
ments is consistent with the mean value within 1.6 stan- 
dard deviations (see Figure 5). We thus conclude that 
GRB 2108124 passes the count ratio test for lensing. 


3.4. x? test 


A simple and robust test that doesn't assume any 
pulse shape was introduced by Nemiroff et al. (2001). 
'This test was recently applied to the claim of Kalantari 
et al. (2021) on GRB 090717 by Mukherjee & Nemiroff 
(2021b). Mukherjee & Nemiroff (2021b) conclude that 
based on the x? test, the claim of gravitational lensing 
can be excluded at the 5o level. 

This test considers the binned lightcurves of the two 
pulses as representing two distributions and asks if 
they are consistent with coming from the same parent- 
distribution. After appropriately re-scaling the first 
pulse and taking into account the background, we per- 
form a x?-test for the hypothesis that the two lightcurves 
are drawn from the same distribution. 

The test statistic is defined the following way: 


a D (FP, (t) — Po(t + At))? 
X £4 FP, (i) + Pat + At) + Bilt) + Bolt + At 
(3) 
where the t is the time, 7 < 1 is the scaling factor, 
Pj; marks the background subtracted counts in the two 
pulses, and B; is the background counts (i = {1, 2}). 
As an example (see Figure 3, top left), we consider 
data from T0 — 2 s to TO +7 s interval (first pulse) and 
compare it with the interval shifted by At — 33.3 s (sec- 
ond pulse). We use 0.256 s resolution, summed over the 
detectors with good signal in the 45-300 keV range, and 
added the signal of the BGO detectors (120-400 keV). 


We use the tte instead of ctime data, because the be- 
ginning of the first pulse has uneven temporal binning 
in the ctime data. 

First, we find the minimum of the x? expression in 
Equation 3 as a function of 7 and we get f = 0.231 (see 
Figure 3, top-left). Next, we calculate the minimum 
x? = 24.3 value for 34 degrees of freedom, which corre- 
sponds to a p-value of 0.89. Thus there is no statistically 
significant difference between the two distributions. 

We compare the lightcurves of the two pulses for dif- 
ferent temporal resolutions, energy ranges and instru- 
ments, in the panels of Figure 3. We consistently find 
there is no statistically significant difference between the 
two pulses. For example, Mukherjee & Nemiroff (2021c) 
performed a preliminary x? analysis on GRB 210812A 
and concluded there is a ~ 2.80 discrepancy between 
the two pulses. Using the same energy range, and detec- 
tor selection as Mukherjee & Nemiroff (2021c) (assum- 
ing in their notation detectors are numbered 1 through 
12, and ctime energy channels 1 through 8) we per- 
formed the x? test on 512 ms resolution lightcurves and 
in the energy range 22-800 keV (Nal detectors only, 
Figure 3, top-right). We find ? = 0.231 minimizes 
x? at a value of yo = 11.5 for 16 degrees of free- 
dom (p-value of 0.78), indicating the two lightcurves 
are consistent with being drawn from the same distri- 
bution. The BGO lightcurve (512 ms resolution) yields 
x? = 14.027(dof = 16) and p-value of 0.597, and for 
ACS (0.4 s resolution) y? = 24.662(dof = 22) and p- 
value of 0.314 (bottom two panels of Figure 3). 


3.5. Bayesian Model Comparison 


Applying an idea from gravitational wave model selec- 
tion Paynter et al. (2021) introduced the Bayesian evi- 
dence to compare the lensing scenario to the case where 
there is no lensing. 

In the no-lens scenario we fit the lightcurve with 
two pulses, with all the parameters left to vary. 
We derive the Bayesian evidence Zn by integrat- 
ing the likelihood over the multi-dimensional param- 
eter space. E.g. for the N1 pulse, this involves 
10 parameters: I(t) = Ini(t|A1, 07,1, 04,1, tmax,1,V1) + 
Ini (t| A2, Or,2, 4,2; tmax,2; V2). In the lensing scenario 
the second pulse is constrained. It has the same 
shape parameters as the first pulse, only differing 
in the normalization (r) and the shift in the peak 
time (At), resulting in 5 + 2 parameters: I(t) = 
Ini (tlA, Or, Oa, ax, V) HIN (4| Afr, 0r, Ca, tmax + At, v). 
Zz, is the evidence in this case. 

Formally, the Bayesian evidence (or simply evidence) 
has the following meaning (e.g. Speagle 2020): from 
Bayes’ Rule, the probability of the model parameters 
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Figure 6. Ratio of counts in the first and second pulse as 
a function of energy and different instruments. The average 
count ratio is 3.45 + 1.03 (horizontal lines.) 


(in our case the peak times, pulse widths, amplitudes, 
etc., denoted by ©), given the observations (D) and a 
model (M, e.g. the lensing scenario using the N1 pulse 
shape) is: P(O|D, M) = P(DJ|O, M)P(O|M)/P(D|M). 
Here P(D|O, M) is the likelihood of the data given the 
model and its parameters and P(O|M) represents the 
our prior knowledge of the parameters. The evidence 
is the denominator in the expression of Bayes’ Rule: 
Z = P(D|M) = Ne P(D|O, M)P(O|M)dO. The inte- 
gral is performed over the hyper-volume Vo constructed 
from all the parameters. 

Calculating the evidence Z is computationally inten- 
sive. It requires integrating the likelihood over a multi- 
dimensional parameter space. We use the bilby python 
package (Ashton et al. 2019) and dynesty nested sam- 
pler (Speagle 2020) to carry out the integration over the 
parameters to find Z. As in Paynter et al. (2021) the 
difference in In Z is the natural logarithm of the Bayes 
Factor (In BF = In Z; — In Zyz) and it can be used to 
decide between the models. The Bayes Factor is addi- 
tive, values from different, independent (e.g. different 
energy ranges) measurements can be added to perform 
model comparison. 

We present the results of this analysis in Table 2. Con- 
veniently, the nested sampling yields also the best fitting 
flux ratio and time delay. Data from every instrument 
and energy range where the second pulse was detectable 
provided positive Bayesian evidence in favor of the lens- 
ing scenario. 


4. GRAVITATIONAL LENS MODELING 


4.1. Point mass lens 
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Figure 7. 1c confidence region of the time delays and the 
flux ratios using the N2 pulse model and bilby nested sam- 
pling across energy ranges and instruments. Gray contour 
lines mark the associated (1 + zi) Mi lens mass values. Note 
the region encompassing At z 33.3 s and r = 4 to 5 is con- 
sistently part of the solutions. 


The simplest mass model is the point mass lens, when 
the mass of the lens is concentrated in a projected region 
smaller than the Einstein radius of the source. In this 
scenario, we can derive the lens mass from the flux ratio 
and the time delay (e.g. Mao 1992): 


3 —1 
CAt (r—1 
1 Mi = — +] 4 
Otam = S (E+) o. 0 
where z is the lens redshift and M is the lens mass, c 
is the speed of light, G is the gravitational constant. 


4.1.1. MCMC lightcurve fitting 


We fit the summed Nal and BGO lightcurves with 
the Markov-chain Monte Carlo method using the emcee 
(Foreman-Mackey et al. 2013) python package with both 
the N1 and N2 pulses (see Figure 1). This method can- 
not select between the lensing and no-lensing scenarios; 
however, it is fast and robust compared to the more 
computation-intensive nested sampling (see section 3.5). 
We can take the result of the lightcurve fits in the lens- 
ing scenario and derive the lens mass, Mi(1 + zz) in the 
point-mass approximation. 

The N1 pulse model leads to a flux ratio of r = 
4.47* 109 and delay time of At; = 33.16*015 s. The 
corresponding point lens mass value is: 


(1+ 2) M; = 1.07* 0618 x 109 Mo. (5) 


For the N2 pulse the flux ratio is r — 410757 and 
the time delay At; = 33.11 + 0.06 s. These lead to a 
point mass lens of mass: 


(1+ zj) M; = 1.13 + 0.06 x 109 Mọ . (6) 
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detector | energy range pulse In(BF) r At 
(keV) model (s) 
Nal 8-45 N1 0.96  414*229 33.8522 
Nal 45-95 N1 QU. 5370. 33.37139% 
Nal 95-300 N1 1.79 5.591797 33.26 + 0.24 
Nal 300-500 N1 0.28 696750 33687333 
BGO 120-400 N1 229. 80771 3316102 
ACS >70 N1 2.89 4.261035 937702 
Swift 25-350 N1 0.50 40225 - 391255 
Nal 8-45 N2 1.94 3.611323 33.25197 
Nal 45-95 N2 2.83 OGG 33.1119: 
Nal 95-300 N2 Bis 535572. 33.301923 
Nal 300-500 N2 Gel 5.51423 33.1155, 
BGO 120-400 N2 5.25 4.531150  33.13+0.23 
ACS >70 N2 5.18 4.204595 33.34 + 0.24 
Swift 25-350 N2 0.55 S72 328 10:05 


Table 2. Bayes factor compilation for different instruments, energy ranges and pulse models using the bilby nested sampling 
method. The flux ratio and the time delay resulting from the fits are shown in the last two columns. 


Nal (8-44 keV) 
Nal (44-95 keV) 
Nal (95-300 keV) 
BGO 

ACS 

BAT 


0.50 0.75 1.00 1.25 1.50 1.75 2.00 
M,(1-zj) (109M,) 


Figure 8. Mass posterior density distribution from the bilby modeling assuming point mass lens. Different colors show different 
energy ranges or instruments indicated in the legend. Continuous lines represent the N2 pulse model, dashed curves represent 
the N1 model. 
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4.2. Singular Isothermal Sphere (SIS) lens model 


This lens model is characterized by the line of sight 
velocity dispersion, oy of its mass distribution. While 
in the point mass lens case, we could constrain the lens 
mass, here, it is only possible to restrict the velocity 
dispersion up to a distance scale D, 


1/4 —1/4 
s -e( 1 cAt Lu) E) kms”. 


32r? D(1+z)r—1 


where D = DorDrs/Dos and Di; mark the angu- 
lar diameter distance combinations between the observer 
(O), lens (L) and source (S), respectively. We assumed 
a GRB redshift z, = 1 and a lens redshift z; = 0.4. 

A simple mass estimate based on the 
virial theorem yields a mas M œ~ 8 x 
10° (c,/15 km s~')?(R/10 pc) Mo, where R = 10 pc 
is an assumed size considered typical for e.g. globular 
clusters for which the SIS model is a good approxima- 
tion. The mass is broadly consistent with the point 
mass lens approximation. 


5. DISCUSSION 


In the previous section, we performed tests to confirm 
GRB 2108124 is affected by strong gravitational lensing. 
Here we discuss the strengths of each test, analyze the 
unlensed GRB properties and present future detection 
prospects. 


5.1. Spectrum 


The most basic spectral test for a lensing scenario 
is the flux ratio test. Because gravitational lensing is 
achromatic, the flux ratio of the two pulses has to re- 
main constant across energy ranges and different instru- 
ments. This is a robust measure because it doesn't de- 
pend on the assumed spectrum. The only caveat to 
consider if the GBM detectors’ pointing has changed 
between the two pulses. In that case, the spectral re- 
sponses change significantly between the pulses and the 
recorded counts cannot be compared across the emission 
episodes of GRB 2108124. The pointing of the detectors 
however has not changed by more than 5 degrees be- 
tween the pulses, which means the detectors’ response 
in the direction of GRB 2108124 is essentially the same 
throughout the duration of the burst. The flux ratios 
across the energy ranges and instruments are clearly 
consistent with being equal (Figure 6). The weighted 
mean is 3.45, and we measure the most significant de- 
viation for the ACS data point (green), which is 1.6 
standard deviations away. 

Mukherjee & Nemiroff (2021c) showed a preliminary 
analysis of the hardness ratio (HR) for ctime channels 


3 and 4 (4 and 5 in their notation) and claimed a 2.2 
c discrepancy between the HR of the two pulses. For- 
mally, the count ratio (CR) test is equivalent to the HR 
test: the hardness ratio of Pulse 1 is HR(P1)=C(P1, 
Ch2)/C(P1, Ch1), where C(P1, Ch2) denotes the count 
rate in pulse 1 for channel 2 (higher energy channel for 
the hardness), and C(P1, Ch1), HR(P2) can be calcu- 
lated analogously. The count ratio in the energy chan- 
nel 1 is CR(Ch1) = C(P1,Ch1)/C(P2,Ch1). Thus from 
CR(Ch1)=CR(Ch2), it follows that HR(P1)=HR(P2). 
Because there are no strong outliers in the count ratio 
test, we expect the data to reflect this in the HR test. 
For the first pulse we find: HR= 1.405 + 0.063 and for 
the second pulse: HR= 1.127 + 0.155. We confirm the 
finding of Mukherjee & Nemiroff (2021c) that the sec- 
ond pulse indeed shows a lower HR. However, taking 
their difference and adding the errors in quadrature, we 
find that the discrepancy is only 1.660, which does not 
invalidate the lensing scenario. 

Next, we carried out a spectral analysis of the two 
pulses using the GBM data. In the Swift-BAT data, the 
second pulse was only visible in the summed lightcurve, 
and INTEGRAL SPI-ACS had only data in 1 energy 
channel. Therefore we only used Fermi-GBM data for 
the detailed spectral analysis. The spectral parameters 
are consistent within errors (Table 1), and a plot of the 
spectral shapes (Figure 4) also shows that the two spec- 
tra overlap when considering the confidence regions. We 
thus conclude that the spectra of the pulses are consis- 
tent with the lensing interpretation. We note the ad- 
vantage of the continuous 128 energy channel data of 
GBM over the 4 channel tte data available for BATSE, 
for which precise spectral fits were not feasible (Paynter 
et al. 2021). 


5.2. Time history 


Gamma-ray instruments have better temporal than 
spectral resolution. E.g., the number of spectral resolu- 
tion elements in 10-1000 keV is S10, while the 5 s dura- 
tion pulse with ~0.1 s resolution has 50 temporal resolu- 
tion elements, where 0.1 s is the approximate timescale 
of variations for GRB 210812A (Bhat et al. 2012). For 
this reason, the temporal study of the lensing scenario 
can provide more constraints. 

The weaker second pulse is closer to the noise level of 
the detectors. Thus the shorter duration of the second 
pulse is in line with expectations for a weaker bursts 
with similar pulse shape. Nonetheless, the duration of 
the two pulses is still consistent within errors as expected 
from a lensing scenario. 

Independent of any pulse models, we first performed 
a x? test to compare the two lightcurves. The x? test 
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determines if the two lightcurves are consistent with 
being drawn from the same distribution. The x? test 
showed that there is no significant difference between 
the lightcurves in different energy ranges and across in- 
struments (Figure 3). We note, however that weakness 
of the second pulse results in relatively large Poisson er- 
rors, and fine temporal structures in the second pulse, if 
there are any, are washed out. This somewhat reduces 
the power of this test, showing only that the general 
shapes of the pulses are consistent. 

Next, we introduced pulse models from the literature 
and fit the lightcurve in different energy bands and dif- 
ferent instruments. Using nested sampling, we evaluated 
the evidence in favor of the lensing scenario using the 
bilby code. Independent of the detector, energy range, 
or pulse model (Table 2), the evidence is consistently for 
the lensing scenario as opposed to the non-lensing sce- 
nario (BF>0 in all cases). We note that the logarithm 
of the Bayes Factor is not necessarily positive for the 
model with less parameters (the lensing model in our 
case). Indeed, e.g. Wang et al. (2021) (their Table 1) 
shows some cases with negative In(BF). 

The Bayes Factor differs depending on which pulse 
model we apply. We find that in energy ranges where the 
second pulse is relatively weak, the evidence for lensing 
is not as strong (but it still favors the lensing). The 
evidence using the N2 pulse model is more compelling 
than in the case of the N1 shape. This can be due to 
the larger number of parameters in the case of the N1 
pulse shape. 

We consider the NaI data alone, and the N2 pulse 
shape fits. The sum of In (Bayes Factors) for the Nal de- 
tectors yields In(BF) ze 10.7. Following Kass & Raftery 
(1995) and Thrane & Talbot (2019) we can assign collo- 
quial meaning to this number. A difference of more than 
8 is considered strong evidence in favor of the lensing 
model. We conclude that the Bayesian evidence thus 
supports the lensing interpretation. We note that the 
BGO and ACS lightcurves also provide additional ev- 
idence and, to a lesser extent, the Swift-BAT data as 
well. 

'The nested sampling provides both a selection crite- 
rion between the models through the Bayesian evidence 
and, at the same time, provides the parameters for time 
delay and flux ratio. We show the values in Table 2 
and Figure 7. Different instruments, detectors, and en- 
ergy ranges all yield a consistent solution, pointing to 
the lensing origin. Most importantly, we can derive the 
lens mass for both pulse shapes and arrive at a consis- 
tent picture indicating a ~ 109 Mc lens (Figure 8). The 
MCMC method yields smaller errors than the nested 


samplig. This is due to the different fitting approaches 
of the two methods (see e.g. Speagle 2020, for details). 


5.9. GRB properties 


If GRB 2108124 had not been lensed, its fluence would 
have been Fo = Fi(1—1) ~ 6.2 x 1078 erg cm™?, 
where Pi is the fluence of the first pulse (see Section 
3.2) and we took r — 4.5 for the numerical value. Using 
this flux value, we can get a broad range for the pos- 
sible redshift of this GRB using empirical correlations 
between gamma-ray properties. 

We scan the 0.1 to 5 redshift range and find that 
z 2, 0.5 is consistent to within one sigma with the Am- 
ati relation (Amati et al. 2002) between the isotropic 
equivalent energy, Eis, and the redshift-corrected peak 
energy, (1 + z)Epeax. While this is a very crude esti- 
mate, it is in line with the average measured redshift of 
GRBs (Bagoly et al. 2006; Jakobsson et al. 2006) and it 
is consistent with our fiducial value of z, = 1. 

The lag-luminosity relationship (Norris et al. 2000; 
Gehrels et al. 2006) provides another estimate of the 
redshift (L ox (Tiag/(1 + 2))~°™). Taking the median 
value of the lag (221.6 ms) and scanning the 0.1 to 5 red- 
shift range, we find the redshift of GRB 210812A within 
the range 0.9 < z < 1.5 is consistent at lo level with 
the lag-luminosity relation. While this is similarly an 
empirical relation, it further reinforces that a redshift of 
Zs = lis a reasonable approximation. 


5.4. Future events 


Even though GRB 210812A had no multiwavelength 
follow-up, as more lensed GRB candidates are observed, 
we expect to have a well-localized counterpart eventu- 
ally. Identifying the afterglow of a lensed GRB, show- 
ing two consistently fading images, would be the smok- 
ing gun evidence for the lensing scenario. The angular 
separation of the two lensed images on the sky are on 
the order of the Einstein radius, which for a 109 M 
point mass is 0E = ((4GM/c)(Drs/DorDos)) ? x 
3 mas(M/106M)!/2(D/0.6 Gpc)~!/? (assuming z; = 
0.4 and zs = 1, mas = ” milli-arcsecond” ). This falls just 
short of the resolution of 10m class optical telescopes 
(x40 mas). The only conceivable way of resolving the 
two images is through very large baseline interferometry 
(VLBI) radio imaging (Casadio et al. 2021). The sen- 
sitivity of VLBI can be as low as Z 10 wJy (Venturi 
et al. 2020) and radio afterglows are detected in signif- 
icant numbers at or above this flux (Chandra & Frail 
2012). VLBI imaging of the two sources will provide 
additional information on the source and lens redshift 
and help constrain the lens model. Capturing a lensed 
GRB with a well-understood lens will allow fully ex- 
ploit the accurate, sub-second time delay measurement 
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achievable for GRBs and will allow for time-delay cos- 
mography. 


5.5. Lensing object 


It is difficult to identify the type of lens that produced 
the two pulses of GRB 210812A. Possible objects include 
black holes or globular clusters. Populations of objects 
can be ruled out based on their number density and 
their contribution to the total lensing probability (Ne- 
miroff 1989). The total number of millilensed GRBs can 
be estimated based on a search for lensing candidates in 
the entire Fermi-GBM GRB catalog. Our goal in this 
paper was to report only on GRB 210812A, and we leave 
population-level studies for a future work. Nonethe- 
less, we can make a few general observations, based 
on e.g. the previously mentioned claims by Kalantari 
et al. (2021) and Wang et al. (2021); Yang et al. (2021). 
For 1-3 lens events and the total number of Fermi-GBM 
GRBs, N > 3100, the lensing rate is (3—9) x 1074. This 
is on par with the rate based on BATSE observations by 
Paynter et al. (2021). 

A black hole mass of M ~ 10°Mo lies at the lower 
end of the supermassive black hole population with mea- 
sured masses (e.g. Woo et al. 2010) and at the upper end 
of the intermediate-mass black hole population (Greene 
et al. 2020). Without detailed counterpart observa- 
tions it is unclear in which group, if any, the lens of 
GRB 210812A belongs to. Further lensed events how- 
ever can provide essential constraints on the rates and 
origin of black holes in this mass range. 

Globular clusters are similarly good lens candidates 
and their masses can indeed reach 109 Mc (Baumgardt 
& Hilker 2018). A SIS model is a good approximation 
to the velocity dispersion in a globular cluster. Paynter 
et al. (2021) found, however, that even globular clusters 
with one order of magnitude smaller mass, ~ 10? Mo, 
do not exist in sufficient numbers to produce the 1/2700 
rate of lensed GRBs for BATSE. In our case, we re- 


quire similar lensing probabilities but with ~ 10°Mo 
globular cluster population. 109 M. globular cluster lies 
above the approximately 2x 10? Mc; turnover mass in the 
mass function (Jordán et al. 2007). This means that the 
larger 105 Mc mass cannot compensate in total lensing 
probability the drop in number density. We thus con- 
clude that the globular cluster lens can be tentatively 
ruled out, and a point mass lens e.g. a black hole is 
more likely. For more definitive statements on the na- 
ture of the lens precise localization and high resolution 
observations will be necessary. 


6. CONCLUSION 


In this paper, we presented multiple lines of evidence 
for GRB 2108124 being gravitationally lensed. The two 
peaks in GRB 210812A have consistent spectrum, time 
profile, and spectral evolution. We determined the flux 
ratio and time delay with multiple methods and arrived 
at a consistent picture. The first pulse is approximately 
4.5 times brighter and the delay between the pulses is 
33.3 s. Assuming a point mass lens, this flux ratio and 
delay corresponds to a lens mass of (1+ z;) M; = 109 Mo. 
There are only a few unchallenged claims in the litera- 
ture for lensed GRB lightcurves. GRB 2108124 presents 
the first strong evidence for lensing a long GRB with a 
flux ratio larger than 2. Future events will benefit from 
high resolution radio observations for definitive proof of 
lensing origin and detailed lens modeling. 
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